#include <graphics.h>

//Course Project
//Olekhno N. 2103/2
//16.11.12

#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <math.h>
#include <time.h>

const int NS=10000;                   //number of steps in modeling
const int N=30;                       //number of layers in body

//uses namespace std;

main()
{
  float Fx[N], Fy[N], dFx[N], dFy[N]; //forces
  float x[N], y[N], dx[N], dy[N];     //coordinates
  float kx, ky;                       //material properties
  float tau;                          //step-time of modeling
  int a;                              //distance between layers
  int i,j,k;
  long t1, t2;                        //for good graphics without OpenGL
  int width, height;                  //window parameters
  
  tau=0.1;                            //initialization if user said nothing
  a=10;
  kx=15;                              //actually here k=k/m
  ky=17;                              //where m is a mass of layer
  y[0]=50;
  x[0]=50;
  Fx[0]=150;
  Fy[0]=100;
  for (i=1; i<N-1; i++)
  {
     y[i]=y[i-1]+a;
     Fy[i]=0;
  }
  for (i=1; i<N-1; i++)
  {
     x[i]=x[i-1];
     Fx[i]=0;
  }
  width=1000;
  height=700;
  
  initwindow(width,height, "Oscillators");
  setcolor(RED);
  
  t1=clock();
  for (k=0; k<NS; k++)                  //step of modeling
  {
     for (i=0; i<N-2; i++)              //calculating new coordinates
     {
        dx[i]= 0.5*Fx[i]*tau*tau;
        x[i]=x[i]+dx[i];
        dy[i]= 0.5*Fy[i]*tau*tau;
        y[i]=y[i]+dy[i];
     }
     dx[N-1]=0;
     dy[N-1]=0;
     
     t2=clock();
     if(t2-t1>40)                     //25 frames per sec
     {
       // cleardevice();                //visualization of movement
        setcolor(WHITE);
        bar(0,0,1000,700);
        setcolor(RED);
        for (i=0; i<N-1; i++)
        {
           moveto(x[i],y[i]);           //yes, passing float to int
           lineto(x[i]+200, y[i]);
        }
        t1=t2;
     }
     
     Fx[0]=Fx[0]-kx*(dx[0]-dx[1]);      //calculating new forces
     Fy[0]=Fy[0]-ky*(dy[0]-dy[1]);
     for (i=1; i<N-1; i++)
     {
        Fx[i]= (dx[i]-dx[i-1])*kx-(dx[i]-dx[i+1])*kx;
        Fy[i]= (dy[i]-dy[i-1])*ky-(dy[i]-dy[i+1])*ky;
     }
  }
  
  getch();
  closegraph();
  //system ("pause");
  return 0;                               
}

